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SUMMARY 

Simplified curve fits for the thermodynamic properties of equilibrium 
air have been devised for use in either the ’’time-dependent" or "shock- 
capturing" computational methods. The accuracies of these curve fits 
are substantially improved over the accuracies of previous curve fits 
appearing in NASA CR-2134. For the "time-dependent" method, curve fits 
were developed for p = p(e, p) , a = a(e, p) , and T = T(e, p) , while for 
the "shock-capturing" method, curve fits were developed for h = h(p, p) 
and T = T(p, p) . The ranges of validity for these curve fits are the 

same as the NASA -ARC RGAS program, namely, temperatures up to 25,000 

-7 3 

and densities from 10 to 10 amagats. These approximate curve fits 
may be particularly useful when employed on advanced computers such as 
the Burrough's ILLIAC IV or the CDC STAR since they avoid the cumbersome 
table- lookup feature of the RGAS program. 
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NOTATION 


a = speed of sound 
e = internal energy 
h = enthalpy 
p = pressure 
R = gas constant 
T = temperature 
7 = h/e 
p = density 

Subscript 

o = standard conditions 



INTRODUCTION 


When computing real gas flows using a finite-difference solution 
of the conservative form of the unsteady Navier-Stokes equations, it 
becomes necessary to determine pressure as a function of density ( p) 
and internal energy (e) . This requirement led to the previous study^ 
in which two different approaches were developed for the case of 

2 

equilibrium air. In the first approach, the NASA-Ames RGAS program 
was modified to allow density and internal energy to be the independent 
variables. This approach permits a very accurate determination of the 
thermodynamic properties of air. Unfortunately, the table-lookup 
feature of the RGAS program is too cumbersome to be effectively employed 
on advanced computers such as the Burrough*s ILLIAC IV or the CDC STAR, 

For this reason, and also to reduce computation time on conventional 
serial computers, simpler approximate methods were investigated in the 
second approach. 

In the second approach, simplified curve fits were devised for 

p = p(e, p) , a=a(e, p) , and T = T(e, p) . In addition, a simplified 

curve fit was made for h = h(p, p) . This latter curve fit is required 

3 

in the ’*shock-capturing" method . The ranges of validity for these corre- 
lation formulas were the same as the RGAS subroutine, namely, temperatures 
up to 25,000 and densities from 10 ^ to 10^ amagats. The accuracies 

of these simplified curve fits were much better than the previous curve 
4 

fits of Barnwell , but they did not approach the accuracy of the modified 
RGAS program. For this reason, the present study was undertaken to sub- 
stantially improve the accuracies of the previous curve fits without 
increasing the required computer time. 
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CONSTRUCTION OF CURVE FITS 


The curve fits were constructed using Grabau-type transition func- 

5 6 4 

tions in a manner similar to Lewis and Burgess and Barnwell . A 

transition function of this type can be used to smoothly connect two 

surfaces y) ^ ~ constant, the Grabau-type 

transition function (with an inflection point) becomes 




Z - fj^(x) + _ X ) 


where K is the parameter which determines the rate at which z changes 
from io f^(x) , and is the location of the inflection point 

as shown in Fig. 1. 



X X 

o 

Fig. 1. Grabau-type transition with 
inflection point. 


In the previous study , the curve fits were constructed by joining 
two Grabau-type transition functions with the equation for a perfect gas. 
In the present study, a substantial improvement in accuracy was achieved 
by joining together as many as five Grabau-type transition functions with 


the perfect gas equation as shown in Fig. 2. 
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The coefficients in the equations for y) and were 

determined using a least squares computer program to fit the data from 
the original NASA RGAS program. The selection of the form of the equations 
for f^(x, y) and ^2^^’ largely a trial-and-error process. By in- 

cluding more terms, a better curve fit was achieved. In fact, if a suffi- 
cient number of terms were retained in f^(x, y) and ^2^^’ accuracy 

of these curve fits could be made to approach that of the RGAS program, 
but with little savings in computer time. 
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EQUATIONS OF CURVE FITS 
P " PCe> fi) 

For the correlation of p = p(e, p) , the ratio y = h/e was curve- 
fitted as a function of e and p so that p can be calculated from 


p = pe(Y - 1) 


( 2 ) 


The general form of the equation used for y was 


'y = a 4- a^Y + a^Z -!- a.YZ -I- a_Y^ + -I- a YZ^ -h 

1ZJ4 3 D / o 

+ 1 + exp [(a ^3 4- a^^Y)(z + 


(3) 


where Y = log p/1 . 292) and Z = log^^ (e/ 78408 . 4) . The units for p 
3 2 2 

are kg/m and the units for e are m /sec . The coefficients a^, a^, ... 

a are given in Table A1 , Appendix A, for the entire range of e and p. 
Id 

It should be noted that many of the terms appearing in Eq. (3) are not 
used over the entire range of e and p. 


a = a(e, P) 


An exact expression for the speed of sound in terms of v was derived 


4 

by Barnwell and may be written as 


a = 


e (Y - 1) 


Y + 






By 


logeP/e, j 


1/2 


( 4 ) 


Because of the errors in the approximate expression for y> Eq. (3) , it 
was found that a much better correlation for a = a(e, p) could be 
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obtained from 


+ ( Y - 1 ) 


Y+ 




3 l°8eVp 


+ K, 




1/2 


log^PLJ 


(5) 


where the coefficients K^, K 2 , and were determined using the least- 
squares-best-fit program in conjunction with the NASA RGAS program. The 
coefficients Kj^, K^, and are tabulated in Table Al, Appendix A. 

T = T(e. p) 

In the calculation of T = T(e, p) , the pressure is first found 
using Eq, (2), and then the temperature is found from the equation 

logj^Q(T/ 151 . 78 ) = bj^ + b^Y + b^Z + b^YZ + b^Z^ + b^Y^ + b^Y^Z 

. ,„2 »9 + »io'' * ‘> 11 " 

+ " 1 . e.p[(b,,Y + b,3)(2 . W^^)J 

where Y = log^Q< p/1 . 225) , X = log^Q<p/ 1.0134 X 10^), and Z = X - Y, 

The units for p are newtons/m^, and the units for T are ^K. The coef- 
ficients b- , b_, b.,, are given in Table A2, Appendix A. These coef- 

1 z lo 

ficients were determined in such a manner as to compensate for the errors 
incurred in the initial calculation of pressure using Eq. (2). 

h = h(p, p) 

For the correlation of h = h(p, p) , the ratio 'y = h/e was curve-fitted 
as a function of p and p so that h can be calculated from 


h = (p/p) 


_Y_ 


Y - 1/ 


(7) 
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The general form of the equation used for y 


Y = 


C- + + c^Z + c,YZ 4- 

12 3 4 


+ c^Y + c^Z + CgYZ 
1 + exp[cg(X + c^qY + ■ 


( 8 ) 


where Y = logj^Q(p/l,292) , X = log^Q(p/1.013 X 10^), and Z = X - Y. The 
coefficients c^, c^, .*«» are tabulated in Table A3, Appendix A. 


T = T(p. p) 


The general form of the equation used for the correlation T = T(p, p) 

was 


+ d2Y + d3Z + d^YZ + d^Z^ 

d^f d^Y^dgZ-Kd^YZ+d^/ 

1 + exp[d^^(Z + d^^)] 

where Y = log^^ ( p/1 . 225) , X = log^Q(p/l.0134 X 10^), and Z = X - Y. 

The coefficients d^, d^, •••, d ^2 given in Table A4, Appendix A. 

For the *' time-dependent" method, the three curve fits p = p(e, p) , 
a = a(e, p) , and T = T(e, p) have been placed in a single subroutine 
named TGAS, The calling sequence and FORTRAN IV listing of this subroutine 
appear in Appendix B. For the "shock-capturing" method, the curve fits 
h = h(p, p) and T = T(p, p) have been placed in separate subroutines, 
each named TGAS. These subroutines could be combined into a single sub- 
routine, if desired, or could be used in their present forms in the same 
computer program if one of the subroutines names were changed. The calling 
sequences and FORTRAN IV listings of these subroutines appear in Appendix C 
and Appendix D, respectively. 
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COMPARISONS WITH RGAS PROGRAM 

Comparisons of the curve fits p = p(e, p) , a = a(e, p) , T = T(e, p) , 
h = h(p, p) , and T = T(p, p) with the original RGAS program are shown in 
Figs. 4, 5, 6, 7, and 8. In order to make the comparisons for the first 
three curve fits, the following procedure was used. First, p and p data 
were supplied, which allowed the original RGAS program to compute e. Then, 
this e and the original p were inputed into the TGAS subroutine to obtain 
p, a, and T. Because of this procedure, pressure is plotted as one of 
the independent variables in Figs, 4, 5, and 6. 

In order to assess the relative accuracies of the present curve fits 

with the RGAS program, 500 data points, along constant density lines ranging 
3 -7 

from 10 to 10 amagats, were selected for a comparison. The maximum 
percentage differences between the RGAS and TGAS programs along each con- 
stant density line are tabulated in Table 1, The accuracies of the present 
curve fits are substantially improved over the accuracies of the previous 
curve fits appearing in NASA CR-2134^. The maximum percentage differences 
for the primary variables p = p(e, p) and h = h(p, p) were found to be 
4.77o and 4.6%, respectively. 

A comparison of the relative computer times required for the TGAS 
subroutines and the NASA RGAS programs on the IBM 360-65 computer are 
given in Table 2. The new TGAS subroutine for finding p = p(e, p) , 
a = a(e, p) , and T = T(e, p) is 2.68 times faster than the modified NASA 
RGAS subroutine, as compared with the old TGAS subroutine which was 2.65 
times faster. These comparisons do not include the time spent by the RGAS 
subroutine in reading the tape. The new TGAS subroutine for finding 
h = h(p, p) is 3.80 times faster than the original RGAS program as com- 
pared with the old TGAS subroutine which was 3.88 times faster, again 







Comparison of curve fits for 
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Fig* 7. Comparison of curve fits for h - h(p, p) • 
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Table 1. Maximum percentage differences between RGAS and TGAS programs. 


Density 

ratio 




Curve 

Fit 














ID 

o 

P 

= p(e, p) a 

= a(e, p) T 

= T(e, p) 

h = 

h(p. 

P) 

T = T(p, p) 

10^ 


2 

,2?o 

2.07c 

3.97o 


1.8% 

2.37o 

10^ 


1 

.3 

1.0 

2.5 


1.4 


3.0 

10^ 


1 

.5 

1.3 

3.3 


1.9 


2.6 

10^ 


1 

.8 

1.1 

2.9 


2.2 


2.0 

10-1 


2 

.5 

2.7 

4.4 


2.4 


2.9 

10-2 


2 

.9 

1.4 

3.0 


2.9 


1.9 

lo"^ 


3 

.7 

3,2 

4.1 


3.3 


2.8 

lO'"^ 


4 

.7 

2.9 

5.3 


3.9 


2.2 

10-^ 


3 

.9 

3.2 

3.1 


4.2 


2.5 

10-^ 


4 

.0 

3.7 

4.0 


3.3 


4.3 

lo"^ 


4 

.2 

5.9 

4.4 


4.6 


3.4 

Table 2. 

Comparison of computer times. 









Number of 

Old 

New 





Curve Fit 

Data Points 

TGAS 

TGAS 



RGAS 


P = P(e, 

P) 





3.73 

sec 

(data 

point 0-1) 









includes tape read 

a = a(e. 

P) 


• 5080 

6,04 sec 

5.98 sec 

16.01 

sec 

(data 

points 2-5080) 

T = T(e, 

P) 





19.74 

sec 









2.80 

sec 

(data point 0-1) 
includes tape read 

h = h(p. 

d) 


5096 

2.61 sec 

2.67 sec 

10.12 

sec 

(data points 2-5096) 







12.92 

sec 
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excluding the tape read time. If there are only a few hundred calls made 
to these real gas subroutines, then the TGAS subroutines are substantially 
faster than the RGAS subroutines when the tape read time is included. For 
instance, if there are 500 calls made to find p *= p(e, p) , a = a(e, p) and 
T = T(e, p) , then the new TGAS is 8.95 times faster than the modified RGAS 
subroutine. 

Comparisons of the values obtained at the juncture points of adjacent 
curve fits (see Fig. 2) are shown in Table 3. The maximum deviations 
between the curve fits at the juncture points of the primary variables 
p = p(e, p) and h = h(p, p) are 0.817o and 0.987., respectively. 

The simplified curve fits developed in this study for the thermodynamic 
properties of equilibrium air allow the user to reduce computer time and 
storage while maintaining good accuracy. This is particularly true in the 
"time-dependent” method, since the simplified curve fits could be used 
until near the end of a calculation when the "steady-state” solution is 
approached. Then, the modified RGAS subroutine could be used to give more 
accurate thermodynamic properties for the final steps. Substantial savings 
in computer time may also result in the "shock-capturing" method since an 
iterative procedure involving h = h(p, p) is required for equilibrium 


calculations . 
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Table 3. Comparison of variables at juncture points. 



Curve 

Fit 


Density 

ratio 

Point A 

Point B 

Point 

: C 

Point 

D 

Point 

E 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 

Lower 

Upper 




10^ 

1.400 

1.393 

1.294 

1.293 

1.228 

1.231 








10^ 

1.400 

1.393 

1.287 

1.286 

1.208 

1.210 








10^ 

1.400 

1.393 

1.280 

1.279 

1.187 

1.189 








10° 

1.400 

1.393 

1.273 

1.272 

1.166 

1.168 








10-^ 

1.400 

1.395 

1.300 

1.298 

1.185 

1.185 

1.153 

1.152 



p 

- P(«, 

P) 

10-2 

1.400 

1.395 

1.284 

1.284 

1.173 

1.167 

1.141 

1.138 






10’^ 

1.400 

1.396 

1.270 

1.271 

1.159 

1.150 

1.127 

1.124 






lO'"* 

1.400 

1.396 

1.257 

1.257 

1.144 

1.132 

1.113 

l.llO 






10‘5 

1.400 

1.397 

1.259 

1.257 

1.131 

1.139 

1.089 

1.094 

1.125 

1.125 




10’* 

1.400 

1.397 

1.245 

1.240 

1.125 

1.128 

1.083 

1.086 

1.U5 

1.U5 




lo"^ 

1.400 

1.398 

1.233 

1.223 

1.119 

1.117 

1.076 

1.077 

1.105 

1.105 




lo’ 

443 

438 

364 

359 

314 

314 








10^ 

443 

438 

356 

352 

296 

297 








10^ 

443 

438 

349 

346 

279 

279 








10° 

443 

438 

341 

339 

260 

260 








10*' 

443 

437 

357 

355 

267 

271 

250 

252 



a 

- •(«. 

P) 

lo"^ 

443 

438 

346 

345 

255 

258 

238 

239 






10"^ 

443 

438 

334 

334 

242 

244 

224 

225 






10’"^ 

443 

438 

323 

323 

227 

229 

210 

211 






io“^ 

443 

440 

294 

316 

221 

226 

186 

191 

218 

224 




10'® 

443 

440 

291 

304 

216 

216 

179 

181 

210 

212 




10*^ 

443 

441 

289 

292 

211 

206 

172 

171 

201 

203 




10^ 

5.73 

5.69 

14.87 

14.95 










lo' 

5.73 

5.69 

14.87 

15.03 










lo' 

5.73 

5.69 

14.87 

15.02 










10° 

5.73 

5.69 

14.87 

15.00 










10-' 

5.74 

5.71 

15.42 

15.93 

42.28 

41.98 

109.2 

106.6 



T 

- T(e, 

p) 

10-2 

5.74 

5.69 

15.36 

15.56 

39.46 

39.17 

98.56 

95.77 






10"^ 

5.74 

5.67 

15.30 

15.21 

36.82 

36.55 

88.92 

88.26 






lO-"^ 

5.74 

5.65 

15.25 

14.86 

34.36 

34.10 

80.23 

82.20 






10-5 

3.79 

3.74 

16.50 

16.50 

28.92 

28.97 

51.05 

53.28 

126.3 

126.9 




10'^ 

3.79 

3.74 

16.11 

16.14 

27.34 

27.37 

48.99 

50.00 

119.2 

119.8 




10-2 

3.79 

3.75 

15.72 

15.79 

25.86 

25.85 

47.00 

46.92 

112.5 

112. L 




10^ 

1,400 

1.404 

1.291 

1.294 

1.253 

1.255 








Iq2 

1.400 

1.404 

1.262 

1.284 

1.231 

1,233 








lo' 

1.400 

1.403 

1.272 

1.274 

1.210 

1.211 








10° 

1.400 

1.403 

1.262 

1.263 

1.188 

1.189 








10-' 

1.400 

1.397 

1.313 

1.312 

1.203 

1.204 

1.155 

1.155 



h 

- h(p. 

P) 

10-2 

1.400 

1.397 

1.287 

1.287 

1,191 

1.183 

1.140 

1.137 






10-' 

1.400 

1.397 

1.262 

1.262 

1.171 

1.163 

1.124 

1.121 






10-"* 

1.400 

1.397 

1.236 

1.238 

1.148 

1. 144 

1.107 

1.106 






10’5 

1.400 

1.384 

1.292 

1.301 

1.159 

1.156 

1.104 

1,104 






10-^ 

1.400 

1.387 

1.260 

1.259 

1.142 

1. 138 

1.093 

1.092 






10-2 

1.400 

1.389 

1.228 

1.217 

1.125 

1,121 

1.083 

1.080 






10^ 

5.720 

5.681 

14.96 

14.91 










1q2 

5.720 

5.681 

14.96 

14.96 










lo' 

5.720 

5.681 

14.96 

15.01 










10° 

5.720 

5.681 

14.96 

15.06 










10-' 

5.736 

5.698 

15.42 

15,86 

42.59 

42.09 

109.2 

109.0 



T 

- T(p, 

P) 

10‘2 

5.736 

5.682 

15.37 

15.57 

39.74 

39.45 

98.77 

98.49 






10-' 

5.736 

5.665 

15.33 

15.28 

37.07 

37.01 

89.34 

88.95 






•*4 














10 

5.736 

5.649 

15.28 

15.00 

34.59 

34.73 

80.80 

80.33 






10-5 

3.810 

3.738 

18.51 

18.47 

36.34 

35.98 

83.50 

83.90 






10"^ 

3.810 

3,738 

18.32 

18.33 

35.60 

35.44 

82.45 

82.45 






10-2 

3.810 

3.738 

18.13 

18.20 

34.92 

34.94 

81.35 

80.95 
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APPENDIX A Coefficients for curve fits 
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Table A3. Coefficients for curve fit h = h(p, p) . 


1 

Density 

Range 


Curve 

Range 


*=2 

'3 

'4 

'5 

•=6 


z 

< 0.30 

1.40000 

0 

0 

0 

0 

0 

Y > -0.50 

0.30 

< Z <1.15 

1.42598 

0.000918 

-0.092209 

-0.002226 

0.019772 

-0.036600 

1.15 

< Z < 1.60 

1.64689 

-0.062155 

-0.334994 

0.063612 

-0.038332 

-0.014468 


Z 

> 1.60 

1.48558 

-0.453562 

-0.152096 

0.303350 

-0.459282 

0.448395 


Z 

< 0.30 

1,40000 

0 

0 

0 

0 

0 


0.30 

< Z < 0.98 

1.42176 

-0.000366 

-0.083614 

0.000675 

0.005272 

-0.115853 

-4.50 < Y < 0.50 

0,98 

< Z < 1,38 

1.74436 

-0.035354 

-0.415045^ 

0.061921 

0.018536 

0.043582 


1.38 

< Z < 2.04 

1.49674 

-0.021583 

-0.197008 

0.030886 

-0.157738 

-0.009156 


Z 

> 2.04 

1.10421 

-0.033664 

0.031768 

0.024335 

-0.178802 

-0.017456 


Z 

< 0.398 

1.40000 

0 

0 

0 

0 

0 


0.398 < Z < 0.87 

1.47003 

0.007939 

-0.244205 

-0.025607 

0.872248 

0.049452 

-7 < y < -4.5 

0.87 

< Z < 1.27 

3.18652 

0.137930 

-1.89529 

-0.103490 

-2.14572 

-0.272717 


1,27 

< Z < 1.863 

1,63963 

-0.001004 

-0.303549 

0.016464 

-0.852169 

-0.101237 


Z 

> 1.863 

1.55889 

0.055932 

-0.211764 

-0.023548 

-0.549041 

-0.101758 


Density 

Range 

Curve 

Range 

"7 

^8 

•=9 

o 

o 

'll 


Z < 0.30 

0 

0 

0 

0 

0 


0.30 < Z < 1.15 

-0.077469 

0.043878 

-15.0 

-1.0 

-1.040 

Y > -0.50 

1.15 < Z < 1.60 

0.073421 

-0.002442 

-15.0 

-1.0 

-1.360 


Z > 1.60 

0.220546 

-0.292293 

-10.0 

-1.0 

-1.600 


Z < 0.30 

0 

0 

0 

0 

0 


0.30 < Z < 0.98 

-0.007363 

0.146179 

-20.0 

-1.0 

-0.860 

-4.50 < Y < 0.50 

0.98 < Z < 1.38 

0.044353 

-0.049750 

-20.0 

-1.04 

-1.336 


1.38 < Z < 2.04 

0.123213 

-0.006553 

-10.0 

-1.05 

-1.895 


Z > 2.04 

0.080373 

0.002511 

-15.0 

-1.08 

-2.650 


Z < 0.398 

0 

0 

0 

0 

0 


0.398 < Z < 0.87 

-0.764158 

0.000147 

-20.0 

-1.0 

-0.742 

-7 < Y < -4.5 

0.87 < Z < 1.27 

2.06586 

0.223046 

-15.0 

-1.0 

-1.041 


1.27 < Z < 1.863 

0.503123 

0.043580 

-10.0 

-1.0 

-1.5U 


Z > Z 1.863 

0.276732 

0.046031 

-15.0 

-1.0 

-2.250 


Table A4. Coefficients for curve fit 
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APPENDIX B 

SUBROUTINE TGAS FOR p = p(e, p) , a = a(e, p) , and T = T(e, p) 

The calling statement for this subroutine is 
CALL TGAS (E, RHO, P, A, T) 

with 

2 2 

E = Internal energy, m /sec 
3 

RHO = Density, kg/m 

, 2 

P = Pressure, newtons/m 
A = Speed of sound, m/sec 
T = Temperature, 

The following logic can be employed when the English system of units is 
desired: 

El = E * 0,0929 
RHOl = RHO * 515.4 
CALL TGAS (El, RHOl, PI, Al, Tl) 

P = PI * 0.02088 
A = Al * 3.281 
T = Tl * 1.80 

with 

2 2 

E = Internal energy, ft /sec 

3 

RHO = Density, slugs/ft 

2 

P - Pressure, Ibs/ft 
A = Speed of sound, ft/ sec 
T = Temperature, °R 
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Listing of TGAS for p = p(e, p) , a = a(e, p) , and T = T(e, p) 


SUBROUTINE TGAS(E,RHO,P,A,TI 

Y2=ALOG10IRH0/1.292l 

Z2=ALOG10(E/78408.4) 

IFfY2.GT.-.50» GO TO 11 
IFlY2.GT.-4.50) GO TO 6 
IF<Z2.GT..65) GO TO 1 
GAHM=1.400 
SNOSO=E*.560 
GO TO 18 

1 IFIZ2.GT.1.50) GO TO 2 

GAMM=1.46543+( .007625+. 000292*Y2)*Y2-< .254500+.017244*Y2) *Z2 
A+< .355907+.015422*Y2-. 163235*Z2I*Z2*Z2 
GAME=2. 304* (-.25450-. 017244*Y?+I .7118 14+.030844*Y2-. 489705*Z2 ) *Z 2) 
GAMR=2.304*(. 00762 5+(-.017244+.015422*Z2)*Z2+.000584*Y2) 

Al=-. 000954 
A2=. 171187 
A3=. 004567 
GO TO 17 

2 IF(Z2.GT.2.20) GO TO 3 
GAS1=2.02636+.058493*Y2 
GAS2=.454886+.027433*Y2 
GAS3=.165265+.014275*Y2 
GAS4=. 1 36685+. 010071*Y2 
GAS5=.058493-.027433*Z2 
GAS6=-. 01427 5+.010071*Z2 
GAS7=EXP(0.285*Y2-30.0*Z2+58.41) 

OERE=-30.0 

DERR=0.285 
Al=. 008737 
A2=. 184842 
A3=-. 302441 
GO TO 15 

3 IF(Z2.GT.3.05) GO TO 4 
GAS 1=1. 60804+. 03479 1*Y2 
GAS2=.188906+.010927*Y2 
GAS3=.l 241 17+. 007277* Y2 
GAS4=.069839+.003985*Y2 
GAS5=.034791-.010927*Z2 
GAS6=-.007277+.003985*Z2 
GAS7=EXP(0.21*Y2-30.0*Z2+80.73) 

0ERE=-30.0 

0ERR=0.21 
Al=. 017884 
A2=. 153672 
A3=-. 930224 
GO TO 15 
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4 IF( Z2.GT.3.38) GO TO 5 
GAS 1=1. 25672+. 007073* Y2 
GAS2=.039228-. 000491* Y2 
GAS3=-.721798-.073753*Y2 
GAS4=-. 198942-.021539*Y2 
GAS5=.007073+. 000491*Z2 
GAS6=.073753-.021539*Z2 
GAS7=EXP(0.425*Y2-50.0*Z2+166.7» 

DERE=-50.0 

DERR=0.325 
Al=. 002379 
A2=. 217959 
A3=. 005943 
GO TO 15 

5 GAMM=-84.0327+(-. 831 761+ . 001 153*Y2 )*Y2+C 72 . 2066+ .491914*Y 2 ) * Z2 
A+(-20.3559-. 07061 7*Y2+1. 90979*Z2 l*Z2**2 

GAME=2.304*(72.2066+.491914*Y2+(-40.7118-.141234*Y2+5.72937*Z2) 

A*Z2» 

GAMR=2.304*(-.831761+.002306*Y2+< .491914-.070617*Z2>*Z2» 

Al=. 006572 
A2=. 183396 
A3=-. 135960 
GO TO 17 

6 IF(Z2.GT..65I GO TO 7 
GAMM=1.400 
SNDS0=E*.560 

GO TO 18 

7 IF<Z2.GT.1.54> GO TO 8 
GASl =1.44813+. 001 292*Y2 
GAS2=.073510+.001948*Y2 
GAS3=-.054745+.013705*Y2 
GAS4=-. 055473+ . 021 874* Y2 
GAS 5=. 001 292-. 001 948* Z 2 
GAS6=-.013705+.021874*Z2 
GAS7=eXP<-10.0*( Z2-1.42) » 

DERE=-1.0 

DERR=0.0 
Al=-. 001973 
A2=. 185233 
A3=-. 059952 
GO TO 15 

8 IF( Z2.GT.2.22) GO TO 9 
GAS1=1.73158+.003902*Y2 
GAS2=.272846-.006237*Y2 
GAS3=-.041419-.037475*Y2 
GAS4=.016984-.018038*Y2 
GAS5=.003902+.006237*Z2 
GAS6=.037475-.018038*Z2 

GAS7=EXPl( -10.+3. 0*Y21*{ Z2-. 025*Y2-2.025I ) 

DERE=3.0*Y2-10.0 

DERR=3.0*Z2+12. 15*Y2-20. 325 

Al=-. 013027 

A2=. 074270 

A3=. 012889 

GO TO 15 
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9 IFf Z2.GT.2.90) GO TO 10 
GAS 1= I- 593 50+. 075324* Y2 
GAS2=.1T6186+.026072*Y2 
GAS3*. 200838*. 058536*Y2 
GAS4=-099687+.025287*Y2 
GAS5=.075324-.026072*Z2 
GAS6=-.058536+.025287*Z2 

GAS7=EXP(-10.0*Z2*(5. 0*Z2-13,5)*Y2*27.01 

DERE=5.0*Y2-10.0 

DERR=5.0*Z2-13.5 

A 1=. 004342 

A2=. 212192 

A3=-. 001293 

GO TO 15 

10 GAS1=1.12688-.025957*Y2 
GAS2=-.013602-.013772*Y2 
GAS3=. 127737*. 087942*Y2 
GAS4=.043104+.023547*Y2 
GAS5=-.025957+.013772*Z2 
GAS6=-. 087942*. 023547*Z2 

GAS7=EXP(-20.0*Z2*(4.0*Z2-13.2I*Y2*66.0) 

DERE=-20.*4.0*Y2 

DERR=4.0*Z2-13.2 

Al=. 006348 

A2=. 209716 

A3=-. 006001 

GO TO 15 

11 IF(Z2.GT..65) GO TO 12 
GAMM=1.400 
SN0SQ=E*.560 

GO TO 18 

12 IF(Z2.GT. 1.681 GO TO 13 
GASl=l.45510-.000102*Y2 
GAS2=.081537-.000166*Y2 
GAS3=-. 128647*. 049454*Y2 
GAS4=-. 101 036*. 03351 8*Y2 
GAS5=-.000102*.000166*Z2 
GA S6=-. 049454*. 03351 8* Z 2 
GAS7=EXP(-15.*(Z2-1.420) 1 
DERE=-15. 

0ERR=0. 

Al=. 000450 
A2=. 203892 
A3=. 101797 
GO TO 15 
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13 IF| Z2.GT.2.A6) GO TO 14 
GAS1=1.59608-.042426*Y2 
GAS2=.192840-.029353*Y2 
GAS3=.019430-.005954*Y2 
GAS4=.026097-.006164*Y2 
GAS5=-. 042426+.029353*Z2 
GAS6= ,005954-. 006164+Z2 
GAS7=EXP(-15.*( Z2-2.050) ) 

DERE=-15. 

OERR=0.0 
Al=-. 006609 
A2=. 127637 
A3=. 297037 
GO TO 15 

14 GAS1=1. 54363-. 049071*Y2 
GAS 2=.l 53562-. 0292 09* Y2 
GAS3=.324907+.077599*Y2 
GAS4=.142408+. 022071*Y2 
GAS5=-.049071+.029209*Z2 
GAS6=-.077599+.02?071*Z2 
GAS7=EXP(-10.0*( Z2-2. 708»» 

DERE=-10.0 

DERR=0.0 
Al=-. 000081 
A2=. 226601 
A3=. 170922 

15 GAS10=l./( 1.+GAS7) 

16 GAS8=GAS3-GAS4*Z2 
GAS9=GAS8*GAS7*GAS10**2 
GAMM=GAS1-GAS2*Z 2-GAS 8*GAS10 
GAME=2.304*(-GAS2+GAS4*GAS10+GAS9*DERE» 

GAMR=2.304*t GAS5<-GAS6*GAS10^GAS9*DERR) 

17 SNDSQ=E* ( A1+ ( GAMM-1 , ) ♦ ( GAMM+A2*GAME I +A3*GAMR ) 

18 A=SORT( SNDSQ) 

P=RH0*E*(GAMM-1. ) 

X2=AL0G10( P/1.0134E+05) 

Y2= Y2+. 0231264 
Z3=X2-Y2 

IF( Y2.GT.-.50) GO TO 29 
Y2.GT.-4.50) GO TO 24 
IF<Z3.GT..30) GO TO 19 
T=P/1287.*RH0) 

RETURN 

19 IF(Z3.GT.1.00) GO TO 20 

T=10**( .2718+.00074*Y2+( .990136-. 004947*Y21*Z3+( .990717 
A+.175194*Y2-( .982407+. 159233* Y2I*Z3» /( l.+EXP« -20.*i Z3-0. 88J ) ) ) 
GO TO 32 

20 IF«Z3.GT.1.35) GO TO 21 

T=10**( 1. 39925+. 167780* Y2+(-. 143168-. 159234*Y2l*Z3+(-.027614 
A-.090761*Y2+( . 307036+ . 1 2 162 1*Y2 I *Z3 I / ( 1 .+EXP< -20.* ( Z3-1. 17) ) ) ) 
GO TO 32 

21 1F(Z3.GT.1.79) GO TO 22 

T = 10**( 1.1 1401+.002221*Y2+( .35l875 + .017246*Y2)*Z3+(-l . 15 099 
A-.173555*Y2+( . 673342+. 088399*Y2) *Z3) /( 1. +EXP ( -20.* { Z3-1 . 56 ) ) ) ) 
GO TO 32 
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22 IFfZ3.GT.2.^7) GO TO 23 

T=10**( 1.01722-.017918*Y2+(.473523«-.025^56*Y2)*Z3+(-2.179 78 
A-.334716*Y2+( .898619+,127386*Y2)*Z3)/( l.+EXP ( -20.*( Z3-2. 22) )) ) 

GO TO 32 

23 T=10**(-45.0871-9.00504*Y2+( 35. 8685+6. 7922*Y2)*Z3-< 6. 77699 
A+ 1. 2737*Y2)*Z3*Z3+ (-.064705+. 025325*Z3)*Y2*Y2) 

GO TO 32 

24 IF(Z3.GT..48) GO TO 25 
T=P/(287.*RH0) 

RETURN 

25 TF(Z3.GT..9165 ) GO TO 26 

T=10**( .284312+. 98791 2*Z3+.001644*Y2) 

GO TO 32 

26 IF(Z3,GT. 1.478) GO TO 27 

T=10**( .502071-.01299*Y2+(. 774818+.025397*Y2)*Z3+{ . 009912 
A-. 150527*Y2+(-.000385+.105734*Y2)*Z3)/(l.+EXP(-15.*( Z3-1.28) ) ) ) 

GO TO 32 

27 IFCZ3.GT.2.176) GO TO 28 

T=10*4| 1. 02294+. 021 53 5*Y2+ (. 4272 13+. 006900* Y2)*Z3+(-. 427 823 
A-.211991*Y2+( .257096+.101192*Y2)*Z3) /( 1. +E XP ( -1 2 . * ( Z3-1 . 778 ) ) ) ) 

GO TO 32 

28 T=10**( 1.47540+. 12962*Y2+(. 254154-.046411*Y2)*Z3+(-.221229 
A-.057077*Y2+( .158116+.03043*Y2)*Z3) /( 1 . +EXP ( 5. *Y2*( Z 3-2 . 40) ) ) ) 

GO TO 32 

29 IF(Z3.GT..48) GO TO 30 
T=P/(287.*RH0) 

RETURN 

30 1F(Z3.GT.1.07) GO TO 31 
T=10**( .279268+. 992172*Z3) 

GO TO 32 

31 T=10**(.2332605-.056383*Y2+( 1.1978 3+. 063121* Y2-.165985*Z 3) *Z3+(-. 8 
A14535+.099233*Y2+( .602385-. 067428*Y2-. 098991*Z3)*Z3)/( l.+EXP( ( 
A5.*Y2-20. )♦( Z3-1.78) ) 1 ) 

32 T= T*151. 777778 
RETURN 

END 



29 


APPENDIX C 

SUBROUTINE TGAS FOR h = h(p, p) 

The calling statement for this subroutine TGAS is 
CALL TGAS (P, RHO, H) 

with 

2 

P = Pressure , newtons/m 
2 

RHO = Density, kg/m 

2 2 

H = Enthalpy, m /sec 

The following logic can be employed when the English system of units is 
desired : 

PI = P/0.02088 
RHOl = RHO * 515.4 
CALL TGAS (PI, RHOl, HI) 

H = Hl/0,0929 

with 

2 

P = Pressure, Ibs/ft 

3 

RHO = Density, slugs/ft 
2 2 

H = Enthalpy, ft /sec 
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Listing of TGAS for h = h(p, p) 

SUBROUTINE TGAS( P t RHO,H) 

Y2= AL0G10(RHO/1.292) 

X2= AL0G10(P/1.013E+05) 

Z3= X2-Y2 

IF (Y2 .GT, -. 50 } GO TO 10 
IF CY2 .GT, -4.50) GO TO 5 
IF (Z3 .GT* •398) GO TO 1 
H= (P/RHO)^3.50 
RETURN 

1 IF (Z3 .GT. 0.870) GO TO 2 
GAS1=1.47003+.OOT939*Y2 
GAS2=.244205+. 025607*Y2 
GAS3=-.872248-.049452*Y2 
GAS4=-.764158+.000147*Y2 
GAS5=eXP(-20.00*( Z3-0.742)) 

GO TO 14 

2 IF (Z3 .GT. 1.270) GO TO 3 
GAS1= 3.18652+.137930*Y2 
GAS2= 1. 89529+. 103490*Y2 
GAS3= 2.14572+.272717^Y2 
GAS4= 2.06586+.223046*Y2 
GAS5=EXP(-15.00^( Z3-1.041)) 

GO TO 14 

3 IF (Z3 .GT. 1.863) GO TO 4 
GAS1= 1. 63963-. 00100436*Y2 
GAS2= .303549-. 0164639*Y2 
GAS3= .852169+. 101237^Y2 
GAS4= .503123+. 0435801^Y2 
GAS5=EXP(-10.00*( Z3- 1.544 ) ) 

GO TO 14 

4 GAS1= 1. 55089+. 0559323*Y2 
GAS2= .211764+. 0235478*Y2 
GAS3= .549041+. 101758*Y2 
GAS4= . 276732+. 0460305*Y2 
GAS5=EXP( -15.00* ( Z3-2 .250) ) 

GO TO 14 

5 IF fZ3 .GT. .300) GO TO 6 
H= (P/RHQ)*3.50 

RETURN 

6 IF (Z3 .GT. 0.980) GO TO 7 
GAS1=1. 421 76-. 000366* Y2 

GA S2=. 0836 14-. 000677* Y2 
GAS3=-.005272+.115853*Y2 
GAS4=-.007363+.146179*Y2 
GAS5=EXP( -20.00* C Z3-0. 860)) 

GO TO 14 

7 IF (Z3 .GT. 1.380) GO TO 8 
GAS 1=1. 74436-. 035354*Y2 
GAS2=.415045-.061921*Y2 
GAS3=-.O10536-.O43582*Y2 
GAS4=.0443534-.049750*Y2 
GAS5=EXP(-20.00*( X2-1 . 04*Y2-1.336) ) 

GO TO 14 
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8 IF (Z3 .GT. 2.0401 GO TO 9 
GASl =1. 49674-. 02 1583*Y2 
GAS2=.197008-,030886*Y2 
GAS3=. 1 57738+. 009 158*Y2 
GAS4=.123213-.006553*Y2 
GAS5=EXP(-10.00*( X2-1.05*Y2- 1.8951 ) 

GO TO 14 

9 GAS1=1.10421-.033664*Y2 
GAS2=-.031768-.024335*Y2 
GAS3=.l 78802 +. 01 7455*Y2 
GAS4=.080373+.002511*Y2 
GAS5=EXP(-15.00*( X2-1 .08+Y2-2.650) I 
GO TO 14 

10 IF(Z3.GT..300) GO TO 11 
H= <P/RHO)*3.50 
RETURN 

11 IFIZ3.GT.1.15J GO TO 12 
GAS 1=1. 42598 + . 00091 8* Y2 
GAS2=. 0922 09+. 002 226*Y2 
GAS3=-. 019772 i+. 036600+Y2 
GAS4=-.0774694+.043878*Y2 
GAS5=EXP(-15.00*( Z3-1.040)) 

GO TO 14 

12 IF(Z3.GT. 1.600) GO TO 13 
GAS1=1. 64689-. 0621547+Y2 
GAS2=.334994-.0636120*Y2 
GAS3=. 038332 2+. 01 44677*Y2 
GAS4=. 0 7342 14-. 002441 7* Y2 
GAS5=EXP(-15.00*( Z3-1.360)) 

GO TO 14 

13 GAS1=1. 48558-. 453562*Y2 
GAS2=.152096-.303350*Y2 
GAS3=. 459282-. 44B395*Y2 
GAS4=.220546-.292293*Y2 
GAS5=EXP(-10.00*( Z3-1.600)) 

14 GAS10= l./(l.+GAS5) 

GAMM=GAS 1-GAS 2*Z3-(GAS3-GAS4*Z3)*GAS 10 
H=(P/RH0»*IGAMM/(GAMM-1.) ) 

RETURN 

END 
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APPENDIX D 

SUBROUTINE TGAS FOR T = T(p, p) 

The calling statement for this subroutine TGAS is 
CALL TGAS (P, RHO, T) 

with 

, 2 

P = Pressure, newtons/m 
3 

RHO = Density, kg/m 

o 

T = Temperature, K 

The following logic can be employed when the English system of units is 
desired: 

PI = P/0.02088 
RHOl = RHO * 515.4 
CALL TGAS (PI, RHOl, Tl) 

T = T1 * 1.80 

with 

2 

P = Pressure, Ibs/ft 

3 

RHO = Density, slugs/ft 
T = Temperature, °R 
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Listing of TGAS for T = T(p, p) 


SUBROUTINE TGA S ( P ,RHG , T ) 

Y2=AL0G10(RH0/1.225) 

X2=ALCG10(F/1.013^E+05) 

Z2=X2-Y2 

IF(Y2.GT.-.50) GO TO 28 
IF(Y2.GT.-A,50) GO TO 23 
IF< Z3.GT. .30> GO TO 19 
T=P/(287.*RH0) 

RETURN 

19 IF( Z3.GT.1.C7) GO TO 20 

T=10**( 2.72964+.003725*Y2+( . 93 8 65 1- . 01 192*Y2 ) *Z3+ ( . 682406+ 
A,C89153*Y2-( . 646 54 1+ . 07 0769* Y2 ) 4 23 1 / « 1 . + EXP ( -20 .♦ ( Z3-0 . 82 ) ) » ) 
GO TO 32 

20 IF( Z3.GT.1.57) GO TO 21 

T = 10**J 2.50246-.04 26 2 7*Y2+( 1. 12924+. 0415 1 7* Y2 )*Z3+ 1 1 , 72067+ 
A.268C08*Y2-( 1.25038+. 179711*Y2)*Z3)/( l.+EXP (-20.*1 Z3-1 .33 )) ) ) 
GO TO 32 

21 IF(Z3.GT.2.24) GO TO 22 

T =10** (2. 4453 1-. 047722* Y2+( 1 . 00488+. 034349*Y2 ) *Z3+< 1 . 95 893+ 
A.316244*Y2-( 1. 0 1 200+ . 15 1 561* Y2» * Z3 » / ( 1. +EXP ( -20.* J Z3-1 . 88 ) ) ) ) 
GO TO 32 

22 T=10**( 2.50342+.026625*Y2+( .838860-.009819*Y2)*Z3+(3.58284+ 

A. 533853*Y2-( 1 . 36 1 4 7+ . 1 95436* Y2 I* Z3 » / ( 1 .+EXP ( -20.* ( Z3-2 .47 ) ) ) ) 

GO TO 32 

23 IF(Z3.GT..48) GO TO 24 
T=P/(287.*RH0) 

RETURN 

24 IFC Z3.GT.. 91651 GO TO 25 

T=10**( . 2616 11 +. 990406* Z3+. 00126 7*Y2) 

GO TO 31 

25 IF(Z3.GT. 1.478) GO TO 26 

T = 10**( .4576 43-. 034 27 2*Y2+(.319119+.046471*Y2)*Z3+»-.073233- 
Q.169816*Y2+( . 04326 4+ . 1 1 1 8 54* Y 2 ) *Z3 ) / ( 1 .+EXP I - 15.*< Z3-1 . 2 8 ) ) )) 
GO TO 31 

26 IF(Z3.GT.2.176) GO TO 27 

T=10**( 1.041 72+. 04196 l*Y2+( .412 752-.0 093 29*Y2)*Z 3+<-.434074- 
A. 196 914*Y2+( . 2648 83+ . 10C599* Y2 ) *Z3) /U .+EXP (-15 .♦ ( Z3~l .778 ) ) ) ) 
GO TO 31 

27 T=10**( .418298-.252100*Y2+(.784048+.144576*Y2)*Z3+(-2.00015- 
A.639022*Y2+( .716053+. 206457* Y2 )♦ Z3) / ( l.+EXP (-10.* ( Z3-2.40) ) ) ) 

GO TO 31 

28 IF(Z3.GT..48) GO TO 29 
T=P/(287.*RH0) 

RETURN 

29 IF( Z3.GT..90) GO TO 30 
T=1C**( . 27407+1. 0C082*Z3) 

GC TO 31 

30 T=10**( .235869-.043304*Y2+( 1.17619+.046498*Y2-.143721*Z3)*Z3+ 
A(-1.3767+.160465*Y2+( 1, C898e-.C83489*Y2-.217748*Z3)*Z3) / 

A(1. + EXP(-10.*(Z3-1.78))) ) 

31 T=T*273.2 

32 T=T/1.8 
RETURN 
END 
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